Ground states and dynamics of multi-component 
Bose-Einstein condensates 



Weizhu Bao * 
Department of Computational Science 
National University of Singapore, Singapore 117543 



Abstract 

We study numerically the time-independent vector Gross-Pitaevskii equa- 
tions (VGPEs) for ground states and time-dependent VGPEs with (or with- 
out) an external driven field for dynamics describing a multi-component Bose- 
Einstein condensate (BEC) at zero or very low temperature. In preparation 
for the numerics, we scale the 3d VGPEs, approximately reduce it to lower di- 
mensions, present a normalized gradient flow (NGF) to compute ground states 
of multi-component BEC, prove energy diminishing of the NGF which pro- 
vides a mathematical justification, discretize it by the backward Euler finite 
difference (BEFD) which is monotone in linear and nonlinear cases and pre- 
serves energy diminishing property in linear case. Then we use a time-splitting 
sine-spectral method (TSSP) to discretize the time-dependent VGPEs with an 
external driven field for computing dynamics of multi-component BEC. The 
merit of the TSSP for VGPEs is that it is explicit, unconditionally stable, time 
reversible and time transverse invariant if the VGPEs is, 'good' resolution in 
the semiclassical regime, of spectral order accuracy in space and second order 
accuracy in time, and conserves the total particle number in the discretized 
level. Extensive numerical examples in 3d for ground states and dynamics 
of multi-component BEC are presented to demonstrate the power of the nu- 
merical methods and to discuss the physics of multi-component Bose-Einstein 
condensates. 
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1 Introduction 



Since its realization in dilute atomic gases [2, 12], Bose-Einstein condensation (BEC) 
of alkali atoms and hydrogen has been produced and studied extensively in the 
laboratory [27], and has afforded an intriguing glimpse into the macroscopic quantum 
world. In view of potential applications, such as the generation of bright beams 
of coherent matter waves (atom laser), a central goal has been the formation of 
condensate with the number of atoms as large as possible. It is thus of particular 
interest to study a scenario where this goal is achieved by uniting two (or more) 
independently grown condensates to form one large single condensate. The first 
experiment involving the uniting of multiple-component BEC was performed with 
atoms evaporately cooled in the |F = 2,m/ = 2) and |1,— 1) spin states of 87 Rb 
[39]. Physically speaking, two independently formed condensates are characterized 
by a random relative phase of their macroscopic wave functions. A "fusing" of two 
condensates thus amounts to locking the relative phase in a dissipative process. 
Currently, there are two typical ways to lock the relative phase: (i). An external 
driven field [39]; (ii). An internal atomic Josephson junction [30]. In fact, recent 
experimental advances in exploration of systems of uniting two or more condensates, 
e.g. in a magnetic trap in rubidium [39] and subsequently in an optical trap in 
Sodium [45], have spurred great excitement in the atomic physics community and 
renewed interest in studying the ground states and dynamics of multi-component 
BEC [27, 18, 30, 20]. 

Theoretical treatment of such systems began in the context of superfluid helium 
mixtures [32] and spinpolarized hydrogen [44], and has now been extended to BEC 
in the alkalis [28, 23, 35, 41]. Theoretical predications of properties of uniting two 
or more condensates, e.g. density profile, dynamics of interacting BEC condensates 
[24], motional damping [30] and formation of vortices [31, 33, 36], can now be com- 
pared with experimental data [27, 3] . Needless to say that this dramatic progress on 
the experimental front has stimulated a wave of activity on both the theoretical and 
the numerical front. In fact, the properties of uniting two or more BEC states at 
temperatures T much smaller than the critical condensation temperature T c [34] are 
usually modeled by the vector Gross-Pitaevskii equations (VGPSs) for the macro- 
scopic vector wave function [40, 34] either with an external driven field [27] or with 
an internal atomic Josephson junction [30]. Note that equations very similar to the 
VGPEs also appear in nonlinear optics where indices of refraction, which depend on 
the light intensity, leads to nonlinear terms like those encountered in VGPEs. 

There have been extensive numerical studies of the time-independent Gross- 
Pitaevskii equation (GPE) for ground states and time- dependent GPE for dynam- 
ics of single-component BEC. For computing ground states of BEC, Bao and Du 
[4] presented a normalized gradient flow (NGF), proved energy diminishing and dis- 
cretized it by a backward Euler finite difference (BEFD) method; Bao and Tang [11] 
proposed a method which can be used to compute the ground and excited states 
via directly minimizing the energy functional; Edwards and Burnett [22] introduced 
a Runge-Kutta type method; other methods include an explicit imaginary-time al- 
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gorithm used in [1] and [19]; a directly inversion in the iterated subspace (DIIS) 
used in Schneider et al. [43] and a simple analytical type method proposes by Dodd 
[21]. For numerical solutions of time-dependent GPE for finding dynamics of BEC, 
Bao et al. [6] presented a time-splitting spectral (TSSP) method, Ruprecht et al. 
[42] used the Crank-Nicolson finite difference method, Cerimele et al. [16, 17] pro- 
posed a particle-inspired scheme. Up to now, only a few numerical simulations on 
multi-component BEC [30, 18, 25]. 

In this paper, we take the 3d vector Gross-Pitaevskii equations (VGPEs) with 
an external driven field for multi-component BEC, make it dimensionless, approx- 
imately reduce it to a 2d VGPEs and a Id VGPEs in certain limits, and discuss 
the approximate ground state solution of VGPEs in (very) weak interaction regime. 
Then we present a normalized gradient flow (NGF) to compute ground states of 
multi-component BEC, prove energy diminishing of the NGF which provides a math- 
ematical justification, discretize it by the backward Euler finite difference (BEFD) 
which is monotone in linear and nonlinear cases and preserves energy diminish- 
ing property in linear case. At last, we use a time-splitting sine-spectral method 
(TSSP), which was studied in Bao et al. [7, 8] for the nonlinear Schrodinger equa- 
tion (NLS) in the semiclassical regime and used for GPE of single-component BEC 
[6], damped GPE for collapse and explosion of BEC [5] and Zakharov system for 
plasma physics [9, 10], to discretize the time-dependent VGPEs with an external 
driven field for computing dynamics of multi-component BEC. The merit of the 
TSSP for VGPEs is that it is explicit, unconditionally stable, easy to program, less 
memory requirement, time reversible and time transverse invariant if the VGPEs is, 
'good' resolution in the semiclassical regime, of spectral order accuracy in space and 
second order accuracy in time, and conserves the total particle number in the dis- 
cretized level. Extensive numerical examples in 3d for ground states and dynamics 
of multi-component BEC are presented to demonstrate the power of the numerical 
methods. 

The paper is organized as follows. In section 2 we start out with the 3d VGPEs 
with an external driven field, make it dimensionless, show how to reduce it to lower 
dimensions. In section 3 we give approximate ground state solution in (very) weak 
interaction regime, present a normalized gradient flow (NGF) to compute ground 
states of multi-component BEC, prove energy diminishing of the NGF, discretize 
it by the backward Euler finite difference (BEFD), as well as apply the NGF and 
its BEFD discretization to a nonlinear two-state model for vortex states dynamics 
in BEC. In section 4, we present the time-splitting sine-spectral (TSSP) method 
for the VGPEs with an external driven field. In section 5 numerical tests of the 
VGPEs for ground states and dynamics of multi-component BEC are presented. In 
section 6 a summary is given. Throughout we adopt the standard / 2 -norm of vectors, 
matrices, and || • || as the standard L 2 -norm for functions, as well as the * operator 
which is used in Matlab for two vectors U = (ui, • • • , %) T and V = (vi, • • • , %) T 
as U * V = (uiVi, ■ ■ ■ , u m vm) T ■ 
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2 The vector Gross-Pitaevskii equations (VGPEs) 



At temperatures T much smaller than the critical temperature T c [34], a BEC for 
M components with an external driven field is well described by the macroscopic 
vector wave function \I> = Vl/(x, t) = (^(x, t), ■ ■ ■ , ^m(x, t)) T whose evolution is 
governed by a self-consistent, mean field vector Gross-Pitaevskii equations [26, 40]. 
If the harmonic trap potential is considered, the VGPEs become 

r^vErfx t) ft 2 ~ ~ ~ 

ih )C = -— V 2 *(x, t) + V(x) * *(x, f) + A(*) * *(x, t) + hf(t)B¥(x, t), 
at 2m 

(2.1) 

where x = (x, y, z) T is the spatial coordinate vector, m is the atomic mass, H = 
1.05 x 10~ 34 [Js] is the Planck constant, f(t) is a given real-valued scalar function, 
B = {bji)Yi = i is a given M xM symmetric real matrix, i.e. bji = bij (j, I — 1, • • • , M), 
V(x) = (Vi(x), • • • , V M (x)) T is the harmonic trap potential, i.e. 

^'( x ) = y (^j - ^o,j) 2 + (y - £ 0j ) 2 + ^ (z - z , 3 ) 2 ) , j = 1, • • • , M 

with (xoj,yoj, Zqj) t and w x j, w y j, uj z j the center and trap frequencies in x, y, and z- 
direction, respectively, of the jth (j — 1, • • • , M) component. For the following we as- 
sume (w.r.o.g.) uj xA < ■ ■ u xM < uj yA ■ ■ ■ < uj ZjM . A(*) = (ii(*), • • • , i M (*)) T 
models the interaction, i.e. 



Aj{9) = Uj! |^i | 2 + ■■■ + U jM \i> 



M 



= — W +••• + — V'm , j = l,---,M, 

m m 

with Uji = 4wh ajl and dji = aij the s-wave scattering length between the jth and Ith 
component (positive for repulsive interaction and negative for attractive interaction, 
j, I — 1, • • • , M). It is necessary to ensure that the vector wave function is properly 
normalized. Specifically, we require 

f |^(x, 0)| 2 rfx = V° > 0, j = 1, • • • , M, (2.2) 

where N® is the number of particles of the jth (j = 1 • • • , M) component at time 
t = 0. 

2.1 Dimensionless VGPEs 

In order to scale the VGPEs (2.1), we introduce 



x a 



3/2 



t = w Xt it, x=— , ijjj(5c,t) — —j==ijjj(x,t), j — 1, - ■ ■ ,M, a 



(2.3) 
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where ao is the length of the harmonic oscillator ground state. In fact, here we choose 
l/u Xi i and ao as the dimensionless time and length units, respectively. Plugging (2.3) 



into (2.1), multiplying by 



— to the jth (j — 1 • • • , M) equation, and then 



removing all ~, we obtain the following dimensionless VGPEs in 3d with an external 
driven field 

dvErfx t) 1 

i m = -2V 2 *(x,t)+V(x) * *(x,t)+A(*) * *(x,t)+/(t)£*(x,t), (2.4) 
where f(t) = f (t/u X;1 ) /uj xA , and 



V(x) 


= (V r i(x), 


^m(x)) t , 








V S (x) 


= 2 ( 7 ij ^ ~ 


x o,i) 2 + 7L' 


(y- 


yoj) 2 + (z - 


^o,j) 2 ) , 


7x,j — 


U x,j 


1 




1 




X 0,j = 


a 


a 


— , J = 1, 

a 


••,M, 


A(*) 


= (^i(*), •• 


•, A M (*)) T 


i 







A,-(*) = ^i|Vi| 2 + --- + ^m|^m| 2 , 







UjiNf _ A-Kh 2 ajl N? _ AirajiNf 



a 



j,Z = 1,.--,M, 



B = Gq 1 B G , 



with 



G = diag(VJV?, 



The VGPEs (2.4) conserves the normalization of the vector wave function 

or the total number of particles 



N(G *)=[ ||G *(x,t)|| 2 2 dx = f; / iV°|^-(x,t)| 2 a'x = iV , t > 0. (2.5) 

When there is no external driven field, i.e. / = in (2.4), the VGPEs (2.4) is time 
reversible, time transverse invariant, and conserves the normalization of the wave 
function for each component or the number of particles of each component 



TV, 



(^) = / l<Mx,t)| 2 dx=l, t>0, j = l, 



M 



and the energy 



with 



N° = N® + • • 



JV M> 



/3 = max |/3,-7 1, 

M 



(2.6) 



(2.7) 



iiv^r+^(x)i^i 2 +^E^i^-i 2 i^ 



dx, j = l,---,M. 
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There are two extreme regimes: one is when (3 <C 1 (<^=^ \f3ji\ <C 1 for all 
j, I — 1, • • • , M), then the system (2.4) describes a weakly interacting condensation. 
The other one is when > 1, then (2.4) corresponds to a strongly interacting 
condensation or to the semiclassical regime or the Thomas-Fermi regime. 



2.2 Reduction to lower dimensions 

In the following two cases, the 3d VGPEs (2.4) without external driven field, i.e. 
/ = 0, can approximately be reduced to 2d or even Id. In the case (disk-shaped 
condensation) 

u> x ,j ~ ujyj « u x>1 , u g j > u Xi i 7 x j « 7 WJ - « 1, 7 ZJ - > 1, j = 1, • • • , M, 

(2.8) 

the 3d VGPEs (2.4) can be reduced to 2d VGPEs with x = (x, y) T by assuming that 
the time evolution does not cause excitations along the z-axis since they have large 
energy of approximately Tuv z j compared to excitations along the x and y-axis with 
energies of about fuv Xi i. Thus we may assume that the condensate wave function 
along the z-axis is always well described by the ground state wave function and set 



V = V 2 (x,y,t) *V 3 (z) with V 3 (z) = ^<f> g (x,y,z) * <$>* g (x,y,z) dxdy 



1/2 



(2.9) 



where * 2 (x, y, t) = (il> 2 ,i(x, y, t), ■ ■ ■ , ^ 2 , M (x : y, t)) T ', V 3 (z) = (^ 3) i(*)» • • • , ^ 3M (z)) T 
and & g (x, y, z) = (0 9jl (x, y,z), • • • , (f> g ,M(x, 2/, -2)) T (see detail in (3.5)) is the ground 
state solution of the 3d VGPEs (2.4) by setting / = and g* denotes the conjugate 
of a function g. Plugging (2.9) into (2.4), then * both sides by &l(z), integrating 
with respect to z over (— oo, oo), we get 

i t] = ~V 2 * 2 (x,t) + (V(x) + C) * * 2 (x,t) + A(* 2 ) * * 2 (x,t), (2.10) 



where 



V(x) = (K(x,?/), V M (x,y)) T , 



y) = \ (txj i x - x o,j? + ih (y - yoj) 2 ) , j = i, • • • , M, 



C = (ci, • • • , c M ) T , 

2 1 /■« |# 3> ,(*) 



c j ^^fj_jz-z , j r\^dz+-j_ i 



dz 



dz, j = l,---,M, 



A(*) = (^i(*), - A M (*)) T , 



= E (&i llM*)| 2 |lM*)| 2 dz) |^,,| 2 , j = 1, • • • , M. 

Since this VGPEs is time-transverse invariant, we can replace \I/ 2 — > * 2 * e _lC */ 2 
which drops the constant vector C in the trap potential and obtain the 2d VGPEs 
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with = ^2 and x 



(x, yf 



. gj^C, t) = _ 1 V 2^ (X) t) + y (x) i t) + A ( ^ i ^ 

The observables are not affected by this. 

Similarly in the case (cigar-shaped condensation) 



(2.11) 



W x j « U y j > W X) i, C^J > U x ,i 7a, j « 1, 7„j > 1, 7 2> j > 1, 1 < j < M, 

(2.12) 

the 3d VGPEs (2.4) can be reduced to Id VGPEs with x = x. Similarly to the 2d 
case, we derive the Id VGPEs 



where 



v(x) = (vi(a;), - v M (x)f, 

V j (x) = ± 1 2 x>j (x-x j) 2 , j = l,---,M, 
A(*) = (^i(*), - A M (*)) T , 

= e ha L \Mv>*)\ 2 llMv,*)| 2 dz 

1=1 v JR 

1/2 



(2.13) 



\<t>gA X iVi Z )\ dx 



J 



M. 



In fact, the 3d VGPEs (2.4), 2d VGPEs (2.11) and Id VGPEs (2.13) with an 
external driven field can be written in a unified way 



0*(x,t) 
where 



--V 2 *(x,t)+V d (x) * *(x,f)+A d (*) * *(x,t)+/(f)B*(x,f), x G I 

(2.14) 



V d (x) = (y d)1 (x), - ^,m(x)) T , 

A d (*) = 4*,m(*)) T , 
A dJ (*) = dJ1 |^| 2 + • • • + /3 dJM |^m| 



,M; 



with 



V, 



d,j 



d=l, 
d = 2, 



Pd,ji - 



2^x,j ( X X 0,j) j 

\ ( x - x o,j) 2 + ilj (y - yo,j) 2 ) , 
\ [ih ( x - x o, 3 ) 2 + il 3 (y - yo, 3 ) 2 + il, 3 (z - z ,) 2 ) , d = 3; 

' Pjilmi \^23,j\ 2 |^23,/| 2 dydz, d=l, 
Piif-oolM 2 Mdz, d=2, 
k Pji, d = 3. 
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The VGPEs (2.14) conserves the normalization of the vector wave function 

or the total number of particles 

r M r 

N(G *)= ||Go*(x,t)||jUx = £ / N^ J (x,t)\ 2 dx = N , t >0. 

JR d ~[ JR d J 

(2.15) 

When there is no external driven field, i.e. / = in (2.14), the VGPEs (2.14) is 
time reversible, time transverse invariant, and conserves the normalization of the 
wave function for each component or the number of particles of each component 

Nj(i/>j)= [ l^(x,f)| 2 <k = l, t>0, j = l,---,M (2.16) 

JR d 

and the energy 

M N 

with 



M 



^|v^l 2 + ^-(x)|^| 2 + ^E^ N 2 M 
z z i=i 



dx, j = 1, • • • , M. 



3 Ground state solution 

To find a stationary solution of (2.14) without external driven field, i.e. / = 0, we 
write 



*(x,t) 



-iUt 



* <E»(x), 



(3.1) 



where U = (yU 1; • • • , Hm) is the chemical potential vector of the multi-component 
condensate and $(x) = (0x(x), • • • , 0m(x)) a real- valued vector function indepen- 
dent of time. Inserting (3.1) into (2.14) gives the following equations for (U, <&): 

W * #(x) = -^A*(x) + V d (x) * *(x) + A d (*) * *(x), x G M d , (3.2) 



under the normalization condition 



dx 



M. 



(3.3) 



This is a nonlinear eigenvalue problem under the constraint (3.3) and any eigenvalue 
vector U can be computed from its corresponding eigenfunction vector $ by 



-|V0,(x)| 2 + V^-(x)|0,(x)| 2 + A^(<&)|0,(x) 



dx 



M 



;|V0,(x)| 2 + ^(x)|0,(x)| 2 + 5:^ |0,(X)| 2 |^-(X) 



1=1 



dx 



l r M 

= ^(*) + -| d Efel0Kx)ri0,(x)| 2 rfx, j = l,---,M. (3.4) 



i=i 
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It is easy to see that critical points of the energy functional Ep(<&) under the con- 
straint (3.3) are eigenfunctions of the nonlinear eigenvalue problem (3.2) under the 
constraint (3.3) and versus versa. In fact, (3.2) can be viewed as the Euler-Lagrange 
equations of the energy functional Ep($) under the constraint (3.3). The multi- 
component BEC ground state solution <& ff (x) is found by minimizing the energy 
Ep(<&) under the constraint (3.3), i.e. 

(V) Find (U g = (fi gA , • • • , h 9 ,m) T , ® g = {4> g ,i, 4> 9 ,m) T e u ) such that 



Ep{& 9 ) = min (J 

1 r M 

to* = = Ef>A*g) + 2 J Rd Efe l^-WI 2 l<M x )| 2 1<J<M, (I 

where the set U is defined as 

U= |* | Ep{*) < oo, J Rd |0,(x)| 2 dx = 1, 1 < j < Mj . 

In non-rotating mult i- component BEC, the minimization problem (3.5) has a 
unique real- valued nonnegative ground state solution «& s (x) > for x G lR d [38]. 
When M = 1, i.e. single-component BEC, the minimizer of (3.5) was computed by 
either a normalized gradient flow [4], or directly minimizing the energy functional 
[11], or the imaginary time method [1, 19], etc. Here we extend the normalized 
gradient flow method for computing ground state solution from single-component 
BEC to multi-component BEC. 



3.1 Normalized gradient flow and energy diminishing 

Consider the following continuous normalized gradient flow 

® t = - V d (x) * $ - A d (*) * $ + U*{t) * x G M. d , t > 0, (3.7) 



$(x, 0) = * (x) = (0o,i( x ) ; • • • , 0o,m(x)) t , 



x G 



!>d. 



(3.8) 



where U<&(t) = • • • , /i#,M (*)) with 



\\M;t)\\ 2 J* 

M 



i|V0,(x,t)| 2 + V^(x) |0,(x,t)| 2 



+ \M^t)\ 2 |0,(x,t)p 



/=1 



rfx, j = l,...,M. (3.9) 



In fact, the right hand side of (3.7) is the same as (3.2) if we view U&(t) as a 
Lagrange multiplier for the constraint (3.3). Furthermore, as observed in [4] for 
single-component BEC, the solution of (3.7) also satisfies the following theorem: 
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Theorem 3.1 Suppose V d (x) > for all x e R d , (3ji > (j,l = 1,---,M) and 
\\4>o,j\\ = lfj = 1,---,M). Then the normalized gradient flow (3.7)-(3.8) is nor- 
malization conservation and energy diminishing, i.e. 

WA;t)f= />;(M)<&HI<M 2 = 1, *>0, j = l,---,M, (3.10) 

^(*) = -E^IIW.0l| 2 = -E^ / R<I W;CM)| 2 *>0, (3.11) 

which in turn implies 

£/?(*(•, h)) > Ep(&{; t 2 )), < t! < t 2 < oo. 



Proof: Multiplying the jth (j = 1, • • • , M) equation in (3.7) by 0-,-, integrating over 
M d , integration by parts and notice (3.9), we obtain 



1 d 

2diJRd TJ 



(x, t) dx — J 4>j d t (pj dx 



I 

-! 

Jr< 



-|V0,(x,t)| 2 + l^-(x)#(x,f) + 



= 0, t>0, j = l,---,M. 

This implies the normalization conservation (3.10). 
Next, direct calculation shows 



rfx 

rfx + ^(t)||0,(-,t)|| 2 

(3.12) 



M 



V0, • v(%) + •2\; / ,(x)o / -v. )( + £ # 



^ iV° 7r« 

M M jyO 



1=1 

M 



dx 



3 

M N 



= E — / 

^ iV° Jk< 



V0, • v(^ 3 ) + 2v dd (pc) ( /> j d t (i> j + E /WjI 2 M^ 



rfx 



M atO M 



*° ,-=1 



+E^ELA,l0,f^dx 



dx 
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N° Jm* 



M 



dx 



= 2 J Rd [-\&<t>j + V d j{*)(t>j + AuWfa 



dt4>j o?x 



™ 2N° r 

= J Rd i- 9 ^ ( x - *) + ^ J (*)& ( x > *)] dt(f>j rfx 

M 2N° d r 

= -E^ii^(^)ir+^(o^/ K j^(x,t)i 2 rfx 

||^(-,t)|| 2 , t>0, 

• • , M) are always real and 
^ / |0,(x,t)| 2 dx = O, j = l,---,M 



_ _J.Ha J. / J-MI2 



(3.13) 



j'=i 

since //*,.,•(£) (j = 1 



due to the normalization conservation. Thus, we easily get 

£/?(*(-, *i)) > £/?(*(-, * 2 )), < t! < t 2 < oo 
for the solution of (3.7). 



□ 



Using argument similar to that in [46], we may also get as t — > oo, <fr approaches 
to a steady state solution which is a critical point of the energy. In non-rotating 
multi-component BEC, it has a unique real valued nonnegative ground state solution 
* s (x) > for all x G R d [38]. We choose the initial data _»„(x) > for x G 
M d , e.g. the approximate ground state solution (3.33) in weakly interacting multi- 
component BEC. Under this kind of initial data, the ground state solution and its 
corresponding chemical potential U g can be obtained from the steady state solution 
of the normalized gradient flow (3.7)-(3.8), i.e. 



<&«,(x) = lim _»(x, t), xG 

t — »oo 



(3.14) 



1 r M 

= = M*g) + 2 J Rd E^|0^(x)| 2 |0 SJ (x)| 2 _x, j = 1, • • • , M. (3.15) 



3.2 Projection 

When one wants to evolve the normalized gradient flow (3.7), (3.8) numerically, 
it is natural to consider the following projection (or splitting) scheme which was 
widely used in physical literatures for computing the ground state solution of single- 
component BEC [4] by constructing a time sequence = t < t x < t 2 < ■ ■ ■ < t n < 
• • • with t n — n k and k > time step: 



$ t = -A* - V d (x) * $ - A d (*) * x G R d , t n < t < t n+1 , n > 0, (3.16) 
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^(x,w) = ^(x,Ci)=iirr#Tii' xGR "> n ^ ' ( 3 - 17 ) 

*(x,0) = * (x), x G M. d ; (3.18) 

where *(x,i±) = (0i(x, t±), ■ ■ ■ , M (x, t±)) T = lim^±$(x,t) and ||0 O j|| = 1 ( j = 
1, • • • , M). In fact, the gradient flow with projection (3.16), (3.17) can be viewed as 
applying the steepest decent method to the minimization problem (3.5) by ignoring 
the constraint $ G U and then projecting back to the set U. The gradient flow 
(3.16) can also be viewed as applying an imaginary time (i.e. t — > —it) in (2.14). The 
normalized step (3.17) is equivalent to solving the following ODE system exactly 

*t(x,i) =U*(t,k) * *(x,t), x G K d , *„<t<t„+i, n> 0,(3.19) 
$(x,t+) = $(x,^ +1 ), x G M. d ; (3.20) 

where W^(t, fc) = ((J,&,i(t, k), • • • , //$,m(£, ^)) T with 

li&,j(t,k) = fi 9J (t n+1 ,k) = -^ ln ll0i(->^+i)H 2 , t n <t< t n+1 , j = 1,---,M. 

(3.21) 

Thus the gradient flow with projection can be viewed as a first-order splitting method 
for the following continuous gradient flow with discontinuous coefficients: 

$ t = - V d (x) * $ - A d (*) * $ + U&(t, k) * x G M d , n > 0, (3.22) 

#(x,0) = $ (x), xGl" 1 ; (3.23) 

Let /c — >■ 0, note (3.21) and (3.16), we see that 

M 

+ l0/(x,t)| 2 |^-(x,op 



i|V^-(x,t)| 2 + \/ dj (x)|0,(x,t)p 



dx, j = l,---,M, (3.24) 



which implies that the problem of (3.22), (3.23) collapses to (3.16), (3.17) as k — > 0. 
Furthermore, using the Theorem 2.1 in [4], we get immediately 

Theorem 3.2 Suppose V d (x) > /or a// x G lR d and ||0 O j|| = 1 (j = l,---,M). 
For Pji = (j, I = 1,---,M), the gradient flow with projection (3.16)-(3.18) is 
energy diminishing under any time step k and initial data <I?o, i- e - 

E {$(; t n+1 )) < E {<f>(; t n )) <■< E {<f>{; 0)) = E Q (<f> ), U = 0, 1, 2, • • • . 

(3.25) 
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3.3 Backward Euler finite difference (BEFD) discretization 

In this subsection, we present a fully discretization of the gradient flow with projec- 
tion (3.16), (3.17) by the backward Euler finite difference (BEFD) which was pro- 
posed in [4] for discretizing a normalized gradient flow for single-component BEC. 

For simplicity of notation we shall introduce the method for the case of one 
spatial dimension (d = 1) with homogeneous Dirichlet boundary conditions. Gen- 
eralizations to d > 1 are straightforward for tensor product grids and the results 
remain valid without modifications. For d—1, the problem becomes 

*t = 2*xx - Vi(x) * * - Ai(*) i *, a < x < b, t n < t < t n+1 , n > 0, (3.26) 
4>j{x,t n+1 ) = 0j-(x,t+ +1 ) = f^!- +1 L a<x<b, n>0, j = l,---,M, (3.27) 

WjV^n+l)]] 

$(x, 0) = & (x), a<x<b, (3.28) 
$(a, t) = t) = 0, t > 0; (3.29) 

with 

\\M\ 2 = I (Pl/x) dx = l, j = 1, • • • , M. 

We choose the spatial mesh size h = Ax > with h = (b — a)/N and N an even 
positive integer, and define grid points by 

Xj :=a + jh, j = 0, 1, • • •, N. 

Let <fr™ = ((0i)", (0m)j) T be the numerical approximation of <&(xj,t n ) = 

(4>i(xj,t n ), , • • •, 4>M(xj,t n )) T . Here we use the backward Euler for time discretiza- 
tion and second-order centered finite difference for spatial derivatives for the gradient 
flow (3.16). The detail scheme is: 



<fj* _ <f>™ i 

3 3 _ x 

k 2h 2 



cf>* „ _ 9** 4- cf>* , 



(0/)r 1= / == j = 0,---,N, Z = 1,.-.,M, n = 0,l,..-, (3.30) 



-V^^-A^)***, 

j = l,---,iV-l, 

*o = ^ = 0, 

'hTZ^Mi ... 
*; = *ofo), i = 0,1, •••,7V. 

It is easy to see that the discretizetion BEFD (3.30) is monotone for any time step 
k > when Vi(x) > and f3ji > = 1, • • • , M). Furthermore, similar to the 
proof of Theorem 3.1 in [4], we can prove the BEFD normalized flow (3.30) is energy 
diminishing for any time step k > when V^x) > and (3ji = (j, / = 1, • • • , M). 

Remark 3.1 Extension of the BEFD discretization ( 3.26) for mutli- component BEC 
can be done as those in the Appendix in [4] for single- component BEC in the cases 
when Vd(x) in 2d with radial symmetry or in 3d with spherical symmetry or cylin- 
drical symmetry, as well as in 2d or 3d for central vortex states. 
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3.4 Approximate ground state solution 

For a weakly interacting condensate, i.e. j3 <C 1 (-<=>■ <C 1, j, / = 1, • • • , M), 
we drop the nonlinear terms (i.e. the last term on the right hand side of (3.2)) and 
find the linear vector Schrodinger equations with the harmonic oscillator potentials 

U * #(x) = -^A#(x) + V d (x) * *(x), x = (xi, • • • , x d ) T e R d , (3.31) 
under the normalization condition (3.3). The ground state solution of (3.31) is [37] 



/& = " 1J 2 ' d '\ J = 1,"--,M, (3.32) 
0^.(x) = ^ Xl ' j ' "J^'^ 1 e -(^iA x i-{ x i)o,j?+---+ix d ,j(x d -{x d ) l)>J ) 2 )/2 _ (3 33) 

It can be viewed as an approximate ground state solution of (3.2) in the case of 
a weakly interacting multi-component BEC. This approximate ground state can 
be used as initial data in the normalized gradient flow (3.7), or (3.16) and (3.17), 
or (3.22) for computing the ground state solution of mult i- component BEC when 
Pji ± 0. 

3.5 Application to a two-state model 

The normalized gradient flow method and its BEFD discretization for multi-component 
BEC can be applied to compute coupled basis wavefunctions with lowest energy of 
the nonlinear two-state model used in [14, 15] for studying vortex dynamics in single- 
component BEC with (or without) an external rotation. For the continence of the 
reader, here we briefly review the derivation of the nonlinear two-state model from 
the Gross-Pitaevskii equation (GPE). Consider the dimensionless GPE for BEC in 
2d with radial symmetry [4, 11, 6]: 



i ij>t(r,6,t) = ~ 



1 d ( dip\ d 2 ij 
r dr \ dr j d6 2 



2 

+ r -ij + P\ij\ 2 ij, (3.34) 



under the normalization condition 

roo r2iv 



rco r2 
JO JO 



{r,6,t)\ 2 r drdO = 1, 



where (r, 6) is the polar coordinate, ip(r, 9, t) is the macroscopic wave function for the 
condensate, (3 is a parameter models the interaction. In order to represent the con- 
densate mean-field wavefunction if) by the superposition of a symmetric component 
(p s and a vortex component <f) v e 10 , we take the ansatz 

V>(r, 9, t) = a s( p s (r- n^e'^ + a^(r; n v )e lt ' e~^\ (3.35) 

where a s and a v are the complex amplitudes of the symmetric and vortex compo- 
nents, respectively. The vortex fraction is < n v = \a v \ 2 < 1 and the symmetric 
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fraction is n s = \a s \ 2 = 1 — n v . The S and (p v are real nonnegative functions, and 
are normalized to unity, i.e. 

2tt / \<f) s (r;n v )\ 2 r dr = 1, 2tt / \<j> v (r; n v )\ 2 r dr = 1. (3.36) 

JO JO 

Plugging (3.35) into (2.10), multiplying both sides by 1 and e~ ie , respectively, and 
then integrating over IR 2 , see detail in [14], we get the following nonlinear two-state 
model 



Hs<f>s = 



Id/ d<f> s 
2r dr \ dr 



+ (3 (n s <j) 2 + 2Fn v <f) 2 v ) <f> s , 



(3.37) 



I d ( d(p t 



2r dr \ dr 



d(f> s (r;n v ) 



dt 



0, (f) v (0;n v ) = 0, lim <p s (r;n v ) = lim (f> v (r;n v ) = 0; (3.39) 



r=0 



where the factor F — 1, but can be adjusted in some cases. In order to study vortex 
dynamics in BEC through the two-state model (3.37), (3.38) [14, 15], one needs to 
find the coupled basis wavefunctions s (r; n v ) and (p v (r; n v ) for any given < n v < 1 
by minimizing the energy E((f) s , <j) v ) defined as 



E(<(> s , <j) v ) = n s E s (((> s , <j> v ) + n v E v ((p s , <f> v ), 

2 



(3.40) 



E s ((j) s , 4>v) = 27T 



E v ((p s ,(p v ) = 2tt 



2 



d<t) s (r]n t 



dr 



+ r (f) s (r;n v ) 



+P (n s (f> 2 s (r; n v ) + 2Fn v (f) 2 v (r; n v )^j </) 2 s (r; n v )^ dr, 



2 



d(j) v (r;n v ) 



dr 



r ' 2 + ^)<fi{r\n v ) 



+f3 (2Fn s <f) 2 s {r; n v ) + n v (f) 2 v {r; n v )) <f) 2 v {r; n v ) dr; 



(3.41) 



under the constraint (3.36). The continuous normalized gradient flow for computing 
the above minimizer is: 

d(f> s (r,t;n v ) 1 d 



dt 



YrTr \ T ~df) -J^s- (3{n s (j) 2 s + 2Fn v (j) 2 v j S + fi 8 (t)(f> a , 



4„(r, t;n v ) 1 d f d<p v \ ( r 2 1 



dt 2r dr \ dr 

{r,t;n v ) 



(j) v - (3 (2Fn s (f) 2 s + n v (j)f) (f) v + /jL v (t)(j) v , 



dt 



\2 2r 2 / 

0, <f> v (0, t; n v ) = 0, lim <f> s (r, t; n v ) = lim </> v (r, t; n v ) = 0, 



r=0 



4> s (r, 0; n v ) = (f> s ,o(r) > 0, <f> v (r, 0; n v ) = (f> v ,o(r) > 0, < r < 00; 



(3.42) 

(3.43) 

(3.44) 
(3.45) 
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with 



and 



/*oo /*oo 

2tt ^ \<p sfi (r)\ 2 r dr = I, 2n J q |0„,o(r)| 2 r dr = 1 



/o°°r |0 s (r,t;n,,)| 2 rfr7o 2 



,(r,t;n„) 



+ r 2 2 (r,t;n„) 



+2(3 (n s 2 (r, i; n„) + 2Fn„0 2 (r, t; 2 (r, i; n v ) dr, 



fj, v (t) = 



J™r\Mr,t;n v )\ 2 drJo 2 



d4> v (r,t;n v ) 



dr 



+ ( r2 + ^) <^( r >*;« 



+2/3 (2Fn s 2 (r, t; n„) + n v (f) 2 v (r, t; n v )) 4> 2 v {r 



dr. 



If we choose the initial data (/> Sj o(r) > and 4> v ,o(r) > for < r < oo, e.g. 
0s,o( r ) = ^T72e _r 7/2 and 0„,o( r ) = ^T72e~ r //2 , then the minimizer 4> s ^{r;n v ) and 
4> Vt g(r; n v ) can be obtained from the steady state solution of the normalized gradient 
flow (3.42)-(3.45), i.e. 



(t>s, g {r;n v ) = lim <f) s (r,t;n v ), (j> v , g (r; n v ) = lim <f) v (r, t; n v ), < r < 00; 

t— >oo t— >oo 



/*oo ^ 



d4> s , g (r;n v ) 



dr 



+ r <j) Syg (r;n v ) 



+2(3 (n s (f) 2 s g (r; n v ) + 2Fn„</> 2 g (r; n v )) 4> 2 %g {r; dr, 



Hv, g = 2vr 



OO <r« 

2 



dr 



+2/5 (2Fn s 2 9 (r; n v ) + n v (f) 2 vg (r; n v )^ g (r; n v ) 



dr. 



Choose R > sufficiently large, denote the mesh size h r = (R — 0)/N, grid 
points r-j = j h r , j = 0, 1, • • • , N, and r _i = (j — |) j — 0, 1, • • • , N + 1. Then 



the BEFD discretization of the normalized gradient flow (3.42)-(3.45) reads: 

1 



j 2 j 2 



2 /i 2 . r,i 

■> 2 



(rj + r,_i) (<^ s )*_i + rj_i(^ s )J_3 



-/3 



j 2 



* 1, j = !,-■■, N, 



ife 



2 ^ r . 



5 _L 



2Fn s 



'(0 S )"_ !+(^r +1 

J 2 J T 2 



n\2 
j > 



(</>„)*, j = l,---,iV-l, 



16 



2 



J 2 



(0 s )t , i = 0, (<f> v )o = (<f) v )* N = 0, 

2 JV ^2 



J 2 



j = 0, ■ • ■ , iV + 1, n = 0,l, 



i = 0,---,JV, n = 0,l,---, 



(&)°! =^,o(r,_i), j = l,---,JV+l, (« .,=(^, 

J 2 J 2 2 2 

(60? = &,,ofa)> J = 0, TV. 



(3.46) 



(3.47) 



Remark 3.2 T/ze normalized gradient flow and its BEFD discretization for the two- 
state model in 2d with radial symmetry can be easily extended to the two-state model 
in [15] in 3d with cylindrical symmetry. 

4 TSSP method for dynamics 

In this section we present a time-splitting sine-spectral (TSSP) method for the 
VGPEs (2.14) with (or without) an external driven field for dynamics of multi- 
component BEC. For simplicity of notation we shall introduce the method in one 
space dimension (d — 1). Generalizations to d > 1 are straightforward for tensor 
product grids and the results remain valid without modifications. For d — 1, the 
equations (2.4) with homogeneous Dirichlet boundary conditions become 

d^f(x t) 1 

i K m ! = --* xx (x, t) + V 1 (x) * t) + Ax(*) * t) + f(t)B*(x, t), 
a<x <b, t > 0, (4.1) 
*(a,t) = *(M) = 0, t>0, (4.2) 
&(x,t = 0) = * (x) = ^ ,M(a:)) T , a < x < 6, (4.3) 



with 



|foMI= f \^oAx)\ 2 dx = 1, j = 1, • • • , M. 

J a 



We choose the spatial mesh size /i = Ax > with h = (b — a)/N for iV an even 
positive integer, the time step k = At > and let the grid points and the time step 
be 

Xj := a + j h, t n :— n k, j = 0, 1, • • • , N, n — 0, 1, 2, • • • 
Let *™ = ((V'i)?) (V'm)™) 71 be the approximation of *(£,-,£„) = (i>i(xj, t n ), 
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From time t — t n to t — t n+ i, the VGPEs (4.1) is solved in three splitting steps. 
One solves first 

. d&(x,t) 1 T . N .„ „ N 

' dt 2 ^ X,t) ( } 

for the time step of length k, followed by solving 

for the same time step, and then by solving 

. d*(x,t) 



dt 



= f(t)B*(x,t), (4.6) 



Equations (4.4) will be discretized in space by the sine-spectral method and inte- 
grated in time exactly. For t G [t n ,t n+ i], the ODE system (4.5) leaves \i/)j(x, t)\ 
(j = 1, • • • , M) invariant in t and therefore becomes 

and thus can be integrated exactly. The solution of (4.7) is 

t) = e -*CVi(*)+A 1 (*(*,t I ,)))(f-f„) i ^ t n < t < tn+1 . (4.8) 

For t G [t n ,t n+ i], the ODE system (4.6) can be solve exactly too. In fact, sine B is 
real and symmetric, there exist an orthonormal real matrix P with P T = P^ 1 and 
a diagonal matrix D = diag(di, • • • , c?m) such that 

B = P D P T . 

Thus the matrix B in (2.4) can be diagonalizable too, i.e. 

B — Gq 1 P D P t G . (4.9) 

Let 

$(x,t) = P T G *(x,t), (4.10) 
Plugging (4.10) into (4.6), noting (4.9), we get 

i^^- = f(t)D9(x,t). (4.11) 

The solution of (4.11) is 

t) = e-' D -C m ds t n ), t n < t < Wi. (4.12) 
Substituting (4.12) into (4.10), we obtain the solution of the ODE system (4.6) 

*0M) = G l Pe- tD fL f(s)ds <S>(x,t n ) 

= Go 1 P e - iD $l S{s)ds P T G *(x,t n ), t n <t<t n+1 . (4.13) 
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From time t — t n to t — t n+ i, we combine the splitting steps via the standard 
second-order splitting: 



N-l 



= E e- ik $l A sin 



1=1 \ 

(2) _ , -i^^O+A^*^))^ ■ ^(1) 

' 3 ' 



iV . 



*f = e 



(3) 



(4) 



-i(Vifo)+Ai(*< 3 >))fc/2 ± ^(3) 



*n + l = 1^^^ (*(4) )| sin 



1 < j < iV - 1; (4.14) 



where ^ = ((^i)/, • • • , {"4>m)i) t (I = 1, • • • , N — 1), the sine coefficients of * with 
*o = *at = 0, are defined as 



In — ^U 1 . (i j I n s 



b — a' 



N 



I = N-l. 



(4.15) 



The overall time discretization error comes solely from the splitting, which is sec- 
ond order in k, and the spatial discretization is of spectral (i.e. 'infinite') order of 
accuracy. It is time reversible and time-transverse invariant if the VGPEs (2.4) is, 
i.e. / = 0. Furthermore, for the stability of the TSSP (4.14), we have the following 
lemma, which shows that the total number of particles in the mult-component BEC 
is conserved for any given real- valued external driven field /, and the number of 
particles of each component is conserved when there is no external driven field, i.e. 

/ = o. 

Lemma 4.1 The time-splitting sine-spectral (TSSP) method (4-14) is uncondition- 
ally stable and conserves the total number of particles in the mult- component BEC. 
In fact, for every mesh size h > and time step k > 0, 



M 



X E^° ll(iM B ll?.= 
\ 1=1 



\ 



M 



N-l 



3=1 



= ||Go*o||p = x E 7V / ) = v ^' n=l,2,---. (4.16) 
\ 1=1 

Furthermore, when f = in (2.4), i.e. without external driven field, we have 



my 



N-l 



||(^)°||i» = l, n = l,---, j = l,---,M. 

(4.17) 
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Proof: From (4.16), noting (4.14), Parsaval equality, we get 

||G * n+1 ||f 2 = ||G * (4) || 2 2 = ||G *( 3 )||f 2 = \\Pe- iD ti +lf{ -' s)ds P T G *W\& 
= ||G * (2) ||?2 = ||G * (1) ||, 2 a = ||G * n ||i 2 3, n = 0,l,---. (4.18) 

The conservation (4.16) is obtained from (4.18) by induction. When / = 0, the 
proof of (4.17) follows the line of the analogous result for the linear Schrodinger 
equation by time-splitting Fourier-spectral approximation in [7, 5]. 

Remark 4.1 When the definite integral Ji*" +1 f(s) ds in (4-14) could not be evalu- 
ated analytically for some very complicated function f , it can be evaluated numeri- 
cally using a numerical quadrature, e.g., the Simpson's rule 

[ tn+1 f(s) ds^^- [f(t n ) + 4/(t n + fc/2) + f(t n+1 )} , j = 0, • • • , N, n > 0. 

Jt n 4 

Remark 4.2 When M = 2, 6 n = 6 2 2 = 0, 612 = b 2 ± = 1 and f(t) = Qcos^Ud t) in 
(2.1) [27], thus f{t) = ^-cos(c<j ( i t/u Xt i) in (2.4). After a simple computation, we 
can get explicitly the solution (4- 13) of the ODE system (4-6) which will be used in 
our numerical example 5 in the next section: 

ipi(x,t) = cos(g(t,t n )) ipi(x,t n ) -ism(g(t,t n )) \Jn$/N? i/j 2 (x,t n ), 
^ 2 (x, t) = -i sm(g(t, t n )) s[n[Jn% ^{x, t n ) + cos(<?(i, t n )) M 
where 

rt rt Q 

9(t,t n ) = / f(s)ds= cos(u d s/u xA ) ds 

= — [sin(uj d t/u X}1 ) - sin(uj d t n /u Xj i)] . 

5 Numerical results 

In this section we report the coupled basis wavefunctions with lowest energy of a 
two-state model and ground states of multi-component BEC by using the normalized 
gradient flow method, and dynamics of multi-component BEC by using the time- 
splitting sine-spectral method. Furthermore we also give a physical discussion on 
our numerical results. 

Example 1 Coupled basis wavefunctions with lowest energy of a two-state model, 
i.e. we choose (3 = 100 and F = 0.79 in (3.37), (3.38). We solve this problem on 
[0, 8], i.e. R = 8 with mesh size h r — ^ and time step k = 0.1 by using the BEFD 
discretization (3. 46), (3. 47). Ths initial data in (3.45) is chosen as 

s ,o(r) = 4n e " r2/2 > Mr) = 472^ rV2 > r * °- 
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The steady state solution is reached when — 0™ 



m+l. 



s v n h2 < e = 10" 



Table 1 displays the values of S (O), energies i? s , E v and -E 1 , chemical potentials /j, s , 
ji v . Figure 1 shows the coupled basis wavefunctions 4> s , g (r) and <f> v ,g(f) for different 
vortex fraction < n v < 1, and Figure 2 shows surface plots of the atomic density 

12 " 2 



function 
vortex fraction C 



a s <Ps, g + a v <f>. 
< n,, < 1. 



with a v = 



and a s = y/1 — n v for different 



n v 


n s 


0s(O) 


E s 


E v 


E 


Us 







1 


0.2381 


3.9459 


NA 


3.9459 


5.7598 


NA 


0.1 


0.9 


0.2517 


3.8697 


5.8901 


4.0717 


5.7939 


6.9516 


0.3 


0.7 


0.2875 


3.7370 


5.4928 


4.2637 


5.8864 


6.6513 


0.5 


0.5 


0.3433 


3.6474 


5.1110 


4.3792 


6.0166 


6.4091 


0.7 


0.3 


0.4450 


3.6291 


4.7618 


4.4220 


6.1946 


6.2276 


0.9 


0.1 


0.7113 


3.7653 


4.4637 


4.3938 


6.4023 


6.0951 


1 





NA 


NA 


4.3689 


4.3689 


NA 


6.0514 



Table 1: Numerical results for a two-state model in 2d in Example 1. 



From Table 1 and Figures 1&2, we can see that, when the vortex state fraction 
n v increases from to 1, the central value of the symmetric state <f> s (0) and the total 
energy E increase, chemical potentials of the symmetric state fi s and vortex state 
fi v increases and decreases, respectively; the atomic density function \ip\ 2 changes 
from ground state (cf. Fig. 2a) to vortex state (cf. Fig. 2f). 

Example 2 Ground state of two-component BEC in 3d with dynamically stable 
inter-component interaction, i.e. a n > 0, a 2 2 > and 011022 — cl\ 2 > [30]. We 
choose M = 2, m = 1.44 x 10 -25 [kg], an = a 2 2 = 5.45 [nm], a\i = a 2 i = 5.24 [nm], 

Vx,l = Vy,l = ^x,2 = Uy,2 = 10 X 27T [l/s], LU Zj i = U z , 2 = y/&V x ,l, ^0,1 = ^0,2 = 

3/0,1 = Vo,2 = Zq,i = zo,2 = 0, / = in (2.1). Plugging these parameters into (2.4), 
we get the dimensionless parameters a = 0.3407 x 10~ 5 [m], /3 n = 0.02010177^°, 
P 12 = 0.0193272A" 2 °, f3 21 = 0.0193272^°, (3 22 = 0.02010177A^ 2 °. We compute the 
ground states of this problem in cylindrical coordinate on (r, z) G VL = [0, 8] x [—4, 4] 
with mesh size h = h r = h z = ^ and time step k — 0.1 by using the BEFD 
discretization for different and A^^. Here we report the results for two cases: 

Case I. A 2 ° = iV°; 
Case II. A 2 ° = 

Table 2 displays the central densities 9j i(O,O) 2 , <p gt2 (0,0) 2 , chemical potential 
^9,1) f^g,2 & n d energy Ep for case I with different iVf . Fi gure 3 shows the ground 
state condensate wave functions for case I. Furthermore Table 3 and Figure 4 show 
similar results for case II. 
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7V° 


0^(0, 0) 




^(0,0) 


A*a,2 


E 





0.5496 


2.4130 


0.5496 


2.4130 


2.4130 


100 


0.4747 


2.7664 


0.4747 


2.7664 


2.5994 


500 


0.3548 


3.6406 


0.3548 


3.6406 


3.1161 


1,000 


0.2969 


4.3481 


0.2969 


4.3481 


3.5650 


3,000 


0.2170 


6.0980 


0.2170 


6.0980 


4.7258 


6,000 


0.1765 


7.7461 


0.1765 


7.7461 


5.8504 


10,000 


0.1513 


9.3204 


0.1513 


9.3204 


6.9388 


20,000 


0.1226 


12.0802 


0.1226 


12.0802 


8.8655 



Table 2: Numerical results for the ground states of two-component BEC in 3d in 
example 2 for case I. 





#,i(0,0) 




$.2(0,0) 


t*9,2 


E 


10 


0.5353 


2.4738 


0.5351 


2.4746 


2.4440 


100 


0.4504 


2.9062 


0.4491 


2.9122 


2.6799 


500 


0.3225 


4.0142 


0.3193 


4.0315 


3.3582 


1,000 


0.2679 


4.8777 


0.2637 


4.9029 


3.9218 


3,000 


0.1963 


6.9773 


0.1902 


7.0200 


5.3429 


6,000 


0.1610 


8.9353 


0.1536 


8.9931 


6.6984 


10,000 


0.1392 


10.7975 


0.1309 


10.8692 


8.0013 


20,000 


0.1145 


14.0520 


0.1053 


14.1471 


10.2956 



Table 3: Numerical results for the ground states of two-component BEC in 3d in 
example 2 for case II. 
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From Tables 2&3 and Figures 3&4, we can see that, when the number of particles 
of the two component are the same, i.e. N® = N®, then the ground state density 
functions for the two components are equal to each other, i.e. 4> 2 g l = <p 2 g2 due 
to the same parameters used for the two components; where when N® ^ N 2 , then 
4>^ i ^ 4>l 2- I n the two cases, when the number of particles in the first condensate N® 
increases, the central value of the density functions 0^(0, 0) and 0^(0, 0) decrease, 
but the the total energy E, and the chemical potentials \i g ^ and u Sj2 increase. 

Example 3 Ground state of two-component BEC in 3d with dynamically unstable 
inter-component interaction, i.e. an > 0, 022 > and 011022 ~ O12 

< [27]. We 

choose M = 2, m = 1.44 x 1(T 25 [kg], a 12 = a 21 = 55.31 = 5.53 [nm], a 1± = 
1.03ai2 = 5.6959 [nm], a 22 = 0.97a 12 = 5.3641 [nm], u Xtl = to z>2 = 47 x 2tt [1/s], 

= Wy,l = Vx,2 = Uy,2 = W Z ,l/V&, ^0,1 = ^0,2 = $0,1 = Vo,2 = 0, / = in (2.1). 

Plugging these parameters into (2.4), we get the dimensionless parameters do = 
0.2643 x 10- 5 [m], (3 n = 0.02708165^°, /3 12 = 0.02629286iV 2 , (3 21 = 0.02629286^°, 
P22 = 0.02550407^°. We compute the ground states of this problem in cylindrical 
coordinate on (r, z) £ VI = [0, 16] x [—12, 12] with mesh size h — h r — h z — ^ and 
time step k = 0.1 by using the BEFD discretization for different N® and N 2 . Here 
we report the results for three cases: 

Case I. £ ,i = zo,2 = 0, N° = N® + N° = 1, 000, 000. Varying the ratio between N® 
and iV 2 °; 

Case II. N$ = N° = 500, 000. £ ,i = -Zq,2 + 0. Varying z 0A ; 

Case II. io,i = -*o,2 = 0.15a , ' N° = Jvf + iV 2 ° = 1,000,000. Varying the ratio 
between N° and V 2 °. 

Table 4 displays the central densities 9i i(O,O) 2 , 9>2 (O,O) 2 , chemical potential 
^3,1; f^g, 2 and energy Ep for case I with different N®. Fi gure 5 shows the ground 
state condensate wave functions for case I. Furthermore Table 5 and Figure 6 show 
similar results for case II and Table 6 and Figure 7 for case III. 



TV? 


iV 2 ° 


#,i(0,0) 




^2(0,0) 


Vg,2 


E 


100,000 


900,000 


0.0007 


47.8143 


0.0453 


47.2278 


33.8714 


300,000 


700,000 


0.0025 


47.9650 


0.0514 


47.2456 


34.0028 


500,000 


500,000 


0.0063 


48.0880 


0.0605 


47.2631 


34.1575 


700,000 


300,000 


0.0129 


48.1999 


0.0759 


47.2779 


34.3323 


900,000 


100,000 


0.0270 


48.3077 


0.1082 


47.2873 


34.5266 



Table 4: Numerical results for the ground states of two-component BEC in 3d in 
example 3 for case I. 

From Tables 4,5&6 and Figures 5,6&7, we can have the following observations: 
(i). In case I, the trap potentials for the two components are the same, when the 



23 





<i(0,0) 




^, 2 (o,o) 


N,2 


E 


0.01 


0.0092 


48.0611 


0.0601 


47.2363 


34.1375 


0.1 


0.0307 


47.3400 


0.0513 


46.4531 


33.5336 


0.5 


0.0370 


43.9392 


0.0425 


43.0278 


30.8042 


2.0 


0.0277 


37.3581 


0.0271 


36.4870 


26.2717 


4.0 


0.0001 


36.7085 


0.0000 


35.8441 


26.0203 



Table 5: Numerical results for the ground states of two-component BEC in 3d in 
example 3 for case II. 





iV 2 u 


#,i(0,0) 




# >2 (0,0) 


N,2 


E 


100,000 


900,000 


0.0023 


44.3215 


0.0452 


47.0379 


33.4594 


300,000 


700,000 


0.0129 


45.9058 


0.0502 


46.5788 


33.1356 


500,000 


500,000 


0.0335 


46.8865 


0.0489 


45.9868 


33.1615 


700,000 


300,000 


0.0463 


47.6052 


0.0285 


45.1929 


33.4916 


900,000 


100,000 


0.0443 


48.1503 


0.0072 


43.9508 


34.1429 



Table 6: Numerical results for the ground states of two-component BEC in 3d in 
example 3 for case III. 



fraction of the number of particles in the first component, i.e. N®/N°, increases, the 
energy E increases, and the chemical potentials for the two components, // 9; i and 
fi 9t 2, increases and decreases, respectively. The reason of this is due to an > a>22- 
Furthermore, we observe a crater in the density function of the first component , 
corresponding to a shell in which the second component reside (cf. Fig. 5e&f). 
This confirms the experimental results (cf. Fig. 1 in [27]). (ii). In case II, when 
the distance between the center of the trap potentials for the two components, i.e. 
£o,i — zq,2, increases, the energy E, chemical potentials for the two components, // 9) i 
and /i 9j2 , decrease. Furthermore, the bigger the distance, the more separation in the 
density functions of the two components (cf. Fig.6). (iii). The above observation 
(i) also holds for case III except that the crater in the density function for the first 
component almost disappears (cf. Fig. 7). This is due to the separation of the 
centers of the trap potentials for the two components. 

Example 4 Dynamics of two-component BEC in 3d with dynamically unstable 
inter-component interaction, i.e. an > 0, 022 > and 011022 — CL\2 < [27]. We 
choose M = 2, m = 1.44 x 10~ 25 [kg], a 12 = a 2 \ = 55. 3A = 5.53 [nm], an = 
1.03ai2 = 5.6959 [nm], a 22 = 0.97a 12 = 5.3641 [nm], u z>1 = u z , 2 = 47 x 2tt [l/s], 

<^x,l — ^y,l — ^>x,2 = <^y,2 = W z ,l/ ^0,1 = ^0,2 = Vo,l = Vo,2 = 0, Z ,l = —^0,2 = 

0.15a , /(x, t) = ilcos(uj d t) in (2.1). We start the simulation with the initial data 
chosen as the ground state of (2.4) computed by setting / = and using BEFD 
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discretization. We take u>d = 6.5 x 2ty [1/s], N® = N$ = 500,000, and solve this 
problem on the box [—16, 16] x [—16, 16] x [—8, 8] with mesh sizes h x = h y = 1/4, 
h z = 1/8 and time step k = 0.0002. Figure 8 shows the time evolution of the mean 
of the density functions for the two components, H^ill 2 , HV^H 2 (noticing that the 
number of particles in the two components are iV°||^i|| 2 , A^H^H 2 , respectively) 
for Q = 6.5 x 2ir [1/s], 65 x 2ir [1/s] and 650 x 2n [1/s]. Furthermore Figure 
9 displays the time evolution of the density functions for the two components for 
n = 65 x 2tt [1/s]. 

The general form of time evolution on the number of particles in the two com- 
ponents is similar for different external driven field frequency Q. When Q is small, 
the number of particles in the first component, i.e. iV°||^i|| 2 , increases, attains its 
peak and then decreases; where the pattern for iV" HV^H 2 is opposite (cf. Fig. 8a), 
which is due to the total number of particles in the two components are conserved. 
When Q becomes bigger, the pattern of iV° ||^i|| 2 oscillates for some time period, at- 
tains its absolute peak, and then oscillates again (cf. Fig. 8b). Initially the density 
functions for the two components are well separated (cf. Fig. 9 first row), around 
at time t = 3.4, the number of particles in the first component attains its maximum 
and a bigger condensate (with approximately 52% bigger in term of the number of 
particles for the first component than that initially at time t — 0) is generated (cf. 
Fig. 8b&9). When Q becomes very big, similar pattern of iV°||^i|| 2 is observed. In 
fact, the bigger Q is, the faster oscillation in the pattern of the number of particles 
in the condensates (cf. Fig. 8a,b&c). 

6 Conclusions 

The ground states and dynamics of multi-component Bose-Einstein condensates at 
temperature T much smaller than the critical condensate temperature T c are stud- 
ied numerically by using the time-independent vector Gross-Pitaevskii equations 
(VGPRs) and time-dependent VGPEs with (or without) an external driven field, 
respectively. We started with the 3d VGPEs, scale it to obtain a dimensionless VG- 
PEs, and showed how to approximately reduce it to a 2d VGPEs and a Id VGPEs 
in certain limits. We provided the approximate ground state solution of the VGPEs 
in the (very) weakly interacting condensates. Then, most importantly, we presented 
a normalized gradient flow (NGF) to compute ground states of multi-component 
BEC, proved energy diminishing of the NGF which provides a mathematical jus- 
tification, discretized it by the backward Euler finite difference (BEFD) which is 
monotone in linear and nonlinear cases and preserves energy diminishing property 
in linear well as we used a time-splitting sine-spectral method (TSSP) to 

discretize the time-dependent VGPEs with an external driven field for computing 
dynamics of multi-component BEC. The merit of the TSSP for VGPEs is that it 
is explicit, unconditionally stable, easy to program, less memory requirement, time 
reversible and time transverse invariant if the VGPEs is, 'good' resolution in the 
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semiclassical regime, of spectral order accuracy in space and second order accuracy 
in time, and conserves the total particle number in the discretized level. Extensive 
numerical examples in 3d for ground states and dynamics of multi-component BEC 
are presented to demonstrate the power of the numerical methods. Finally, we want 
to point out that equations very similar to the VGPEs are also encountered in non- 
linear optics. In the future we plan to apply the powerful numerical methods to 
study vortex states and their dynamical stability in mult i- component Bose- Einstein 
condensates. 
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1 2 3 + 5 

a) 




2 + 6 

b) 

Figure 1: The coupled basis wavefunctions for a two-state model in example 1 for 
different vortex fraction n v = 0,0.1,0.3,0.5,0.7,0.9,1 (in the order of decreasing 
peak), a). Symmetric state (f> s , g (r), b). Vortex state (f> v ,g(f)- 
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Figure 2: Surface plot of the atomic density function \ip\ 2 for different vortex fraction 
< n v < 1. a), = 0, b). n v = 0.1, c). n„ = 0.3, d). n v = 0.5, e). n„ = 0.9, f). 
n„ = 1.0. 
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cj. ~ 

Figure 3: Ground state solution in 3d with cylindrical symmetry in ex- 
ample 2 for case I. Condensate wave function on two lines for N® = 
0, 100, 500, 1000, 3000, 6000, 10000, 20000 (in the order of decreasing peak): a). On 
the line z = <f> g> i(r, 0) = (fi g ,2( r , 0); b). On the line r = (f> g ,i(0, z) = 9i 2(O, z). c). 



Surface plot of the condensate density function |0 9i i| 



|0 5j2 | 2 for N° = 20000. 
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_e). 

Figure 4: Ground state solution in 3d with cylindrical symmetry in ex- 
ample 2 for case II. Condensate wave function on two lines for N® = 
10,100,500,1000,3000,6000,10000,20000 (in the order of decreasing peak). On 
the line z = 0: a). fl ,i(r, 0); c). 3j2 (r, 0). On the line r = 0: b). fl ,i(O, z); d). 
0s,2(O,2). Surface plot of the condensate density functions for N® = 20000: e). 



2 ; f). 
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*1 g t'.°) 



0J04 




;p yinn " i;; " :f 5 ' 



-5 



e). - ° o i f). 

Figure 5: Ground state solution in 3d with cylindrical symmetry in example 3 for 
case I. Condensate wave function on two lines for N^/N = 0.1,0.3,0.5,0.7,0.9 (in 
the order of increasing at the origin). On the line z = 0: a). (f> g ,i(r, 0); c).(j) gi 2{r, 0). 
On the line r = 0: b). ff ,i(O, z); d). 9i2 (O,z). Surface plot of the condensate 
density functions for = JV 2 ° = 500,000: e). \<p g ,i\ 2 ; f). |0 9 , 2 | 2 - 
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Figure 6: Ground state solution in 3d with cylindrical symmetry in example 3 
for case II. Condensate wave function on two lines for £0,1/^0 = —-20,2/00 = 
0.01,0.1,0.5,2.0,4.0. On the line z = 0: a). 9 ,i(r,O); c).0 ffi2 (r, 0). On the line 
r = 0: b). c/)g i(0, z); d). (f) g2 (0,z). Surface plot of the condensate density functions 

u 2 - 



for 2o >:l /ao = -^0,2/^0 = 2.0: e). |0 9 ,i| 2 ; f). 
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Figure 7: Ground state solution in 3d with cylindrical symmetry in example 3 for 
case III. Condensate wave function on two lines for N?/N° = 0.1,0.3,0.5,0.7,0.9. 
On the line z = 0: a). fl ,i(r, 0); c).0 3j2 (r, 0). On the line r = 0: b). s ,i(O, z) (in the 
order of decreasing peak); d). <f>g j2 (0,z) (in the order of increasing peak). Surface 
plot of the condensate density functions for N® = = 500,000: e). |0 3 i| 2 ; f). 
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Figure 8: Time evolution of the mean of the density functions for the two compo- 
nents, HV'ill 2 , H^H 2 (labeled as ||-Ei|| 2 and H-E^H 2 ; respectively, which indicates the 
time evolution of the number of particles in the two components, iV°||^i|| 2 , iV^H^H 2 , 
respectively) for different driven field frequencies f2. a). Q — 6.5 x 2n [1/s]; b). 
65 x 2tt [1/s]; c). 650 x 2tt [1/s]. 
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Figure 9: Time evolution of the density functions for the two components for the 
driven field frequencies Q = 65 x 2tt [1/s] at different time, from top to bottom: t = 
0.0,0.24,0.58,0.98,1.42,1.96,2.52,3.4. Left column: |^i| 2 ; middle column: |^ 2 | 2 , 
right column: iV'il 2 + 1^2 1 2 - 
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